#include "rectCells.hpp"
#include "rectMesh.hpp"

namespace hms
{

namespace rectCells_internal
{

auto areas( const RectMesh& mesh ) -> AreasExpression {
	return (
		mesh.dx().matrix() *
		mesh.dy().matrix().transpose()
	).array().reshaped();
}

auto circumfs( const RectMesh& mesh ) -> CircumfsExpression {
	static constexpr Index ng { 1 }; // nGhostLayers
	const Index
		  nxg { mesh.nCellsX()+2*ng },
		  nyg { mesh.nCellsY()+2*ng };

	return 2 * (
		mesh.dx().replicate(1, nyg) +
		mesh.dy().transpose().replicate(nxg, 1)
	).reshaped();
}

auto centroids( const RectMesh& mesh ) -> CentroidsExpression {
	const auto& xCoords { mesh.cells().xCoords().matrix() };
	const auto& yCoords { mesh.cells().yCoords().matrix() };

	return concat_vertical(
		xCoords.transpose().replicate( 1, yCoords.size() ),
		yCoords.transpose().replicate( xCoords.size(), 1 ).reshaped().transpose()
	);
}

} // namespace rectCells_internal


RectCells::RectCells( const RectMesh& mesh ) :
	StructCellsBase<RectCells>(mesh)
{}

auto RectCells::areasImp() const -> AreasReturnType {
	#ifdef RECTMESH_PRECALC_AREAS
		return m_areas;
	#else
		return rectCells_internal::areas( mesh() );
	#endif
}

auto RectCells::circumfsImp () const -> CircumfsReturnType {
	#ifdef RECTMESH_PRECALC_CIRCUMFS
		return m_circumfs;
	#else
		return rectCells_internal::circumfs( mesh() );
	#endif
}

auto RectCells::centroidsImp() const -> CentroidsReturnType {
	#ifdef RECTMESH_PRECALC_CENTROIDS
		return m_centroids;
	#else
		return rectCells_internal::centroids( mesh() );
	#endif
}

void RectCells::storeData(){

	auto centroidCoords = [](const auto& vertexCoords){
		return
			(vertexCoords.head( vertexCoords.size() - 1) +
			 vertexCoords.tail( vertexCoords.size() - 1) ) / 2;
	};
	m_xCoords = centroidCoords( mesh().xCoords() );
	m_yCoords = centroidCoords( mesh().yCoords() );

	auto centroidDiffs = [](const auto& cCoords){
		return 
			cCoords.tail(cCoords.size() - 1) -
			cCoords.head(cCoords.size() - 1);
	};
	m_dx = centroidDiffs(m_xCoords);
	m_dy = centroidDiffs(m_yCoords);

	#ifdef RECTMESH_PRECALC_AREAS
		m_areas = rectCells_internal::areas( mesh() );
	#endif

	#ifdef RECTMESH_PRECALC_CIRCUMFS
		m_circumfs = rectCells_internal::circumfs( mesh() );
	#endif

	#ifdef RECTMESH_PRECALC_CENTROIDS
		m_centroids = rectCells_internal::centroids( mesh() );
	#endif

}
/**
 * @brief Returns the x-coordinates of cell centroids, including ghost cells
 * 
 * @return const ArrayXs& 
 */
auto RectCells::xCoords() const -> const ArrayXs& { return m_xCoords; }
/**
 * @brief Returns the y-coordinates of cell centroids, including ghost cells
 * @see RectCells::xCoords()
 * 
 * @return const ArrayXs& 
 */
auto RectCells::yCoords() const -> const ArrayXs& { return m_yCoords; }

/** @brief Returns the distance in x of all cell centroids including ghost cells
 * as an array. */
auto RectCells::dx() const -> const ArrayXs& { return m_dx; }

/** @brief Returns the distance in y of all cell centroids including ghost cells
 * as an array. */
auto RectCells::dy() const -> const ArrayXs& { return m_dy; }

} // namespace hms
